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(57) Abstract: A method for statistically 
reconstructing an X-ray computed 
tomography image produced by a single 
X-ray CT scan having a polenergetic 
source spectrum and an image 
reconstructor which utilize a convergent 
statistical algorithm which explicitly 
accounts for the polyenergetic source 
spectrum are provided. First and second 
related statistical iterative methods for 
CT reconstruction based on a Poisson 
statistical model are described. Both 
methods are accelerated by the use of 
ordered subsets, which replace sums over 
the angular index of a sinogram with a 
series of sums over angular subsets of the 
sinogram. The first method is generalized 
to model the more realistic case of 
polyenergetic computed tomography 
(CT). The second method eliminates beam 
hardening artifacts seen when filtered 
back projection (FBP) is used without 
post-processing correction. The method 
are superior to FBP reconstruction in 
terms of noise reduction. The methods are 
superior to FBP reconstruction in terms of 
noise reduction. The method and image 
reconstructor of the invention are effective 



in producing corrected images that do not suffer from beam hardening effects. 
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STATISTICALLY RECONSTRUCTING AN X-RAY COMPUTED TOMOGRAPHY IMAGE WITH BEAM 
HARDENING CORRECTION 

BACKGROUND OF THE INVENTION 
5 1. Field of the Invention 

The present invention relates to statistical methods for reconstructing 
a polyenergetic X-ray computed tomography image and image reconstructor 
apparatus and, in particular, to methods and reconstructor apparatus which 
reconstruct such images from a single X-ray CT scan having a polyenergetic source 
10 spectrum. 

2. Background Art 

X-ray computed or computerized tomography (i.e. CT) provides 
structural information about tissue anatomy. Its strength lies in the fact that it can 
provide "slice" images, taken through a three-dimensional volume with enhanced 
15 contrast and reduced structure noise relative to projection radiography. 

Figure 1 illustrates a simple CT system. An X-ray source is 
collimated and its rays are scanned through the plane of interest. The intensity of 
the X-ray photons is diminished by tissue attenuation. A detector measures the 
photon flux that emerges from the object. This procedure is repeated at sufficiently 

20 close angular samples over 180° or 360°. The data from different projections are 
organized with the projection angles on one axis and the projection bins (radial 
distance) on the other. This array is referred to as the sinogram, because the 
sinogram of a single point traces a sinusoidal wave. Reconstruction techniques have 
the goal of estimating the attenuation map of the object that gave rise to the 

25 measured sinogram. 
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Figures 2a-2c illustrate the evolution of CT geometries. Figure 2a 
is a parallel-beam (single ray) arrangement, much like what was found in a first- 
generation CT scanner, The major drawback of this arrangement is long scan time, 
since the source detector arrangement has to be translated and rotated. 

5 The fan-beam geometry of Figure 2b reduces the scan time to a 

fraction of a second by eliminating the need for translation. By translating the 
patient table as the source detector arrangement rotates, one gets an effective helical 
path around the object leading to increased exposure volume and three-dimensional 
imaging. 

10 The latest CT geometry is the cone-beam arrangement, shown in 

Figure 2c. It further reduces scan time by providing three-dimensional information 
in one rotation. It is most efficient in its usage of the X-ray tube, but it suffers from 
high scatter (* 40%). It is also the most challenging in terms of reconstruction 
algorithm implementation. 

15 Two dominant effects, both a function of the X-ray source spectrum, 

govern tissue attenuation. At the lower energies of interest in the diagnostic region, 
the photoelectric effect dominates. At higher energies, Compton scattering is the 
most significant source of tissue attenuation. 

The linear attenuation coefficient jx(x,y,z,E) characterizes the overall 
20 attenuation property of tissue. It depends on the spatial coordinates and the beam 
energy, and has units of inverse distance. For a ray of infinitesimal width, the 
mean photon flux detected along a particular projection line Li is given by: 

E[Y ( ] = f I.{E)e^ h{mE)dl dE (D 

♦ 

where the integral in the exponent is taken over the line L i and I;(E) incorporates the 
25 energy dependence of the incident ray and detector sensitivity. The goal of any CT 
algorithm is to reconstruct the attenuation map \i from the measured data [Y,,... ,Y N ] 
where N is the number of rays measured. 
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Filtered Back Projection 

Filtered back projection (FBP) is the standard reconstruction 
technique for X-ray CT. It is an analytic technique based on the Fourier slice 
theorem. 

5 Use of the FFT in the filtering step of FBP renders the algorithm 

quite fast. Moreover, its properties are well understood. However, because it 
ignores the noise statistics of the data, it results in biased estimators. It also suffers 
from streak artifacts when imaging objects with metallic implants or other high- 
density structures. 

10 Polvenersetic X-rav CT 

In reality, the attenuation coefficient \i is energy dependent and the 
X-ray beam is poly energetic. Lower energy X-rays are preferentially attenuated. 
Figure 7 shows the energy dependence of the attenuation coefficients of water 
(density 1.0 gm/cm 3 ) and bone (at density 1.92 gm/cm 3 ). A hard X-ray beam is one 
15 with higher average energy. Beam hardening is a process whereby the average 
energy of the X-ray beam increases as the beam propagates through a material. 
This increase in average energy is a direct consequence of the energy dependence 
of the attenuation coefficient. 

With a polyenergetic source, the expected detected photon flux along 
20 path Li is given by (1). If one were to ignore the energy dependence of the 
measurements and simply apply FBP to the log processed data, some attenuation 
map m would be reconstructed that is indirectly related to the source spectrum and 
object attenuation properties. 

Beam hardening leads to several disturbing artifacts in image 
25 reconstruction. Figures 8 and 9 show the effect of beam hardening on the line 
integral in bone and water. In the monoenergetic case, the line integral increases 
linearly with thickness. With a polyenergetic beam, the soft tissue line integral 
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departs slightly from the linear behavior. The effect is more pronounced for high 
Z (atomic number) tissue such as bone. 

This non-linear behavior generally leads to a reduction in the 
attenuation coefficient. In bone, beam hardening can cause reductions of up to 
5 10%. Thick bones also generate dark streaks. In soft tissue, the values are 
depressed in a non-uniform manner, leading to what has been termed the "cupping" 
effect. In addition, bone areas can "spill over" into soft tissue, leading to a 
perceived increase in the attenuation coefficient. 

Ream Hardenins Correction Methods 

10 Because of SNR considerations, monoenergetic X-ray scanning is not 

practical. Beam hardening correction methods are therefore necessary for 
reconstructing artifact-free attenuation coefficient images. An ideal reconstruction 
method would be quantitatively accurate and portable to different scanning 
geometries. It would somehow reconstruct ii(x,y,E), retaining the energy 

15 dependence of the attenuation process. This is difficult, if not impossible, to 
achieve with a single source spectrum. A more realistic goal is to remove or reduce 
the beam hardening artifacts by compensating for the energy dependence in the data. 

There are a wide variety of schemes for beam hardening artifact 
reduction. Existing methods fall into three categories: dual-energy imaging, 
20 preprocessing of projection data and post-processing of the reconstructed image. 

Dual-energy imaging has been described as the most theoretically 
elegant approach to eliminate beam hardening artifacts. The approach is based on 
expressing the spectral dependence of the attenuation coefficient as a linear 
combination of two basis function, scaled by constants independent of energy. The 
25 two basis functions are intended to model the photo-electric effect and Compton 
scattering. This technique provides complete energy dependence information for 
CT imaging. An attenuation coefficient image can, in principle, be presented at any 
energy, free from beam hardening artifacts. The method's major drawback is the 
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requirement for two independent energy measurements. This has inhibited its use 
in clinical applications, despite the potential diagnostic benefit of energy 
information. Recently, some work has been presented on the use of multi-energy 
X-ray CT for imaging small animals. For that particular application, the CT 
5 scanner was custom built with an energy-selective detector arrangement. 

Commercial beam hardening correction methods usually involve both 
pre-processing and post-processing, and are often implemented with a parallel or 
fan-beam geometry in mind. They also make the assumption that the object consists 
of soft tissue (water-like) and bone (high Z). Recently, these methods were 
10 generalized to three base materials and cone-beam geometry. 

Pre-processing works well when the object consists of homogeneous 
soft tissue. Artifacts caused by high Z materials such as bone mandate the use of 
post-processing techniques to produce acceptable images. 

The attenuation coefficient of some material k is modeled as the 
15 product of the energy-dependent mass attenuation coefficient m k (E) (cm 2 /g) and the 
energy-independent density p(x,y) (g/cm 3 ) of the tissue. Expressed mathematically, 

M(x,y,E)= £ m k (E)p k (x,y)r k (x,y) (2) 

where K is the number of tissue types in the object and ^(x.y) = 1 if (x,y) e tissue 
k and r*(x,y) = 0 otherwise. 

20 For the classical pre-processing approach, the object is assumed to 

consist of a single tissue type (K= J) and to have energy dependence similar to that 
of water, i.e., 



p(x,y,E) = m^E^iXyy) 

* m w (E)p'°"{x,y) 
-5- 
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where rr^fE) is the mass attenuation coefficient of water and /f* is the effective soft 
tissue density. One can rewrite (1) as follows: 

E[Y,)= J I, (E)e l * " { "' E)a dE 

= f I l {E)e m ' (E) ^ p '' Hz - y) dE 

where 7, is the line integral of the density along path L,. Each function F t (T^ 
i=l, ma .,N is 1-1 and monotone decreasing and hence invertible. The goal of the 
10 pre-processing method is to estimate {7^}^ and from that reconstruct (using FBP) 
an estimate p(x,y) of the energy-independent density //*. In other words, 
p(x,y) = FBP{f f } where f, = F t ' l (Y f ) . This pre-processing approach is inaccurate 
when bone is present, but is often the first step in a post-processing bone correction 
algorithm. 

15 Post-processing techniques first pre-process and reconstruct the data 

for soft tissue correction, as explained above. The resulting effective density image 
is then segmented into bone and soft tissue. This classification enables one to 
estimate the contributions of soft tissue and bone to the line integrals. These 
estimates are used to correct for non-linear effects in the line integrals. The final 

20 artifact-free image is produced using FBP and displays density values independent 
of energy according to the following relationship: 



(5) 
(6) 
(7) 



P(x,y) = P sof '(x,y) + kP boM {*>y) W 



where X 0 is some constant independent of energy that maintains image contrast. 
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Although post-processing accomplishes its goal of eliminating energy 
dependence, it suffers from quantitative inaccuracy. As explainer elsewhere, the 
parameter X 0 is somewhat heuristically estimated. In addition, applying post- 
processing to 3-D or cone-beam geometries can be computationally expensive. 
5 Another beam hardening correction of interest is known. This algorithm is iterative 
(but not statistical). At each pixel, it assumes that the attenuation coefficient is a 
linear combination of the known attenuation coefficients of two base materials, and 
it iteratively determines the volume fractions of the base materials. The algorithm 
depends on an empirically-determined estimate of effective X-ray spectrum of the 
10 scanner. The main limitations of this approach is that the spectrum estimate 
captures the imaging characteristics for a small FOV only and that prior knowledge 
of the base materials at each pixel is necessary. 

Statistical Reconstruction for CT 

Many of the inherent shortcomings of FBP are adequately (and 
15 naturally) compensated for by statistical methods. 

Statistical methods are a subclass of iterative techniques, although the 
two terms are often used interchangeably in the literature. The broader class of 
iterative reconstruction techniques includes non-statistical methods such as the 
Algebraic Reconstruction Technique (ART) which casts the problem as an algebraic 
20 system of equations. Successive substitution methods, such as Joseph and Spital's 
beam-hardening correction algorithm, are also iterative but not statistical. Hence, 
statistical methods are iterative, but the opposite is not necessarily true. 

Statistical techniques have several attractive features. They account 
for the statistics of the data in the reconstruction process, and therefore lead to more 
25 accurate estimates with lower bias and variance. This is especially important in the 
low-SNR case, where deterministic methods suffer from severe bias. Moreover, 
statistical methods can be well suited for arbitrary geometries and situations with 
truncated data. They easily incorporate the system geometry, detector response, 
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object constraints and any prior knowledge. Their main drawback (when compared 
to FBP) is their longer computational times. 

Statistical reconstruction for monoenergetic CT was shown to 
outperform FBP in metal artifact reduction, in limited-angle tomography and to 
5 have lower bias-noise curves. 

All single-scan statistical X-ray reconstruction algorithms assume 
(either implicitly or explicitly) monoenergetic X-ray beams, and thus do not deal 
with the issue of beam hardening artifacts. The future potential of statistical 
methods to correct for beam-hardening was speculated by Lange et al. in 1987, but 
10 no methods have been proposed for realizing this potential. One exception is the 
area of dual-energy CT. Dual-energy systems operate based on the principle that 
the attenuation coefficient can be expressed as a linear combination of two energy 
basis functions and are capable of providing density images independent of energy. 

For X-ray CT images, with typical sizes of 512 x 512 pixels or 
15 larger, statistical methods require very long computational times. This has hindered 
their widespread use. 

The article of Yan et al. in Jan. 2000 IEEE Trans. Med. Im. uses 
a polyenergetic source spectrum, but it is not statistical and it does not have any 
regular izat ion. There is no mathematical evidence that their algorithm will 
20 converge. 

In general, all previous algorithms for regularized statistical image 
reconstruction of X-ray CT images from a single X-ray CT scan have been either: 

1) based explicitly on a monoenergetic source assumption, or 

2) based implicitly on such an assumption in that the 
25 polyenergetic spectrum and resulting beam hardening effects 

were disregarded. 
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SUMMARY OF THE INVENTION 

An object of the present invention is to provide a method for 
reconstructing a polyenergetic X-ray computed tomography image and an image 
reconstructor apparatus, both of which utilize a statistical algorithm which explicidy 
accounts for a polyenergetic source spectrum and resulting beam hardening effects. 

Another object of the present invention is to provide a method for 
reconstructing a polyenergetic X-ray computed tomography image and an image 
reconstructor apparatus, both of which utilize a statistical algorithm which is 
portable to different scanner geometries. 

In carrying out the above objects and other objects of the present 
invention, a method for statistically reconstructing a polyenergetic X-ray computed 
tomography image to obtain a corrected image is provided. The method includes 
providing a computed tomography initial image produced by a single X-ray CT scan 
having a polyenergetic source spectrum. The initial image has components of 
materials which cause beam hardening artifacts. The method also includes 
separating the initial image into different sections to obtain a segmented image and 
calculating a series of intermediate corrected images based on the segmented image 
utilizing a statistical algorithm which accounts for the polyenergetic source spectrum 
and which converges to obtain a final corrected image which has significandy 
reduced beam hardening artifacts. 

The step of calculating may include the steps of calculating a gradient 
of a cost function having an argument and utilizing the gradient to minimize the cost 
function with respect to its argument. The step of calculating the gradient may 
include the step of back projecting. The cost function preferably has a regularizing 
penalty term. 

The step of calculating the gradient may include the step of 
calculating thicknesses of the components. The step of calculating thicknesses may 
include the step of reprojecting the segmented image. 
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The step of calculating the gradient may include the step of 
calculating means of data along paths and gradients based on the thicknesses of the 
components. 

The argument may be density of the materials at each image voxel. 

5 The method may further include calibrating the spectrum of the X-ray 

CT scan. 

The method may further include displaying the final corrected image. 

The step of calculating the gradient may include the step of utilizing 
ordered subsets to accelerate convergence of the algorithm. 

10 Further in carrying out the above objects and other objects of the 

present invention, an image reconstructor apparatus for statistically reconstructing 
a polyenergetic X-ray computed tomography image to obtain a corrected image is 
provided. The apparatus includes means for providing a computed tomography 
initial image produced by a single X-ray CT scan having a polyenergetic source 

15 spectrum wherein the initial image has components of materials which cause beam 
hardening artifacts. The apparatus further includes means for separating the initial 
image into different sections to obtain a segmented image and means for calculating 
a series of intermediate corrected images based on the segmented image utilizing a 
statistical algorithm which accounts for the polyenergetic source spectrum and 

20 which converges to obtain a final corrected image which has significantly reduced 

» 

beam hardening artifacts. 

The means for calculating may include means for calculating a 
gradient of a cost function having an argument and means for utilizing the gradient 
to minimize the cost function with respect to its argument. The cost function 
25 preferably has a regularizing penalty term. 
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The means for calculating the gradient may include means for back 

projecting. 

The means for calculating the gradient may include means for 
calculating thicknesses of the components. 

5 The means for calculating thicknesses may include means for 

reprojecting the segmented image. 

The means for calculating the gradient may include means for 
calculating means of data along paths and gradients based on the thicknesses of the 
components . 

10 The argument may be density of the materials at each image voxel. 

The means for calculating the gradient may include means for 
utilizing ordered subsets to accelerate convergence of the algorithm. 

The above objects and other objects, features, and advantages of the 
present invention are readily apparent from the following detailed description of the 
15 best mode for carrying out the invention when taken in connection with the 
accompanying drawings. 

BRIEF DESCRIPTION OF THE DRAWINGS 
FIGURE 1 is a schematic view of a basic CT system; 

FIGURES 2a-2c are schematic views of various CT geometries 
20 wherein Figure 2a shows a parallel-beam (single ray) arrangement, Figure 2b shows 
a fan-beam geometry and Figure 2c shows a cone-beam arrangement; 

FIGURE 3 is a schematic view which illustrates system matrix 
computation for the fan-beam geometry; 

-11- 



WO 02/067201 



PCT7US01/04894 

FIGURE 4 shows graphs which illustrate convex penalty functions; 



FIGURE 5 shows graphs which illustrate the optimization transfer 

principle; 

FIGURE 6 shows graphs which are quadratic approximations to the 
5 Poisson log likelihood; 

FIGURE 7 shows graphs which illustrate attenuation coefficient 
energy dependence; 

FIGURE 8 shows graphs which illustrate beam hardening induced 
deviation of line integral from linearity in water; 

10 FIGURE 9 shows graphs which illustrate beam hardening induced 

deviation of line integral from linearity in bone; and 

FIGURE 10 is a block diagram flow chart of the method of the 
present invention. 

BEST MODE FOR CARRYING OUT THE INVENTION 

15 In general, the method and system of the present invention utilize a 

statistical approach to CT reconstruction and, in particular, iterative algorithms for 
transmission X-ray CT. 

Although the method and system of the invention are described herein 
for a single-slice fan-beam geometry reconstruction, the method and system may 
20 also be used with cone-beam geometries and helical scanning. The method and 
system may also be used with flat-panel detectors. 
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Statistical Reconstruction for X-ray CT 

The physical and statistical models for the problem of X-ray CT 
reconstruction is described herein, and an objective function is obtained. 
Maximizing the objective function by some appropriate iterative algorithm yields 
5 the reconstructed image. 

Monoenerpetic Model and Assumptions 

For the purposes of describing the basic principles underlying this 
invention, for the benefit of skilled practitioners who wish to implement the method, 
the simplified setting in which the transmission source is monoenergetic is first 
10 reviewed. Later the description is generalized to the more realistic polyenergetic 
case. 

The image in object space (attenuation coefficient) is parameterized 
using square pixels. The goal of the algorithm becomes to estimate the value of the 
discretized attenuation coefficient at those pixels. Let \i = .,|ip]' be the vector 
15 of unknown attenuation coefficients having units of inverse length. 

■ 

The measurements in a photon-limited counting process are 
reasonably modeled as independently distributed Poisson random variables. In 
transmission tomography, the mean number of detected photons is related 
exponentially to the projections (line integrals) of the attenuation map. The 
20 measurements are also contaminated by extra background counts, caused primarily 
by scatter in X-ray CT. Thus, the following model is assumed for measurements: 

Y, - Poisson{b i e' [A ' tJ i4-r i },i=l,...,N (9) 

where = Ii(E 0 ) and N is the number of measurements (or, equivalently, the 
number of detector bins). The notation [Ap] { = eJL, a^ ^ represents the ith line 
25 integral. The N x p matrix A = {aij} is the system matrix which accounts for the 
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system geometry as well as any other significant physical effects such as detector 
response. 

Figure 3 illustrates one method for computing the elements of A in 
the fan-beam case. For ray i and pixel j f a (J is the normalized area of overlap 
5 between the ray beam and the pixel. The term r { is the mean number of background 
events, b { is the blank scan factor and Y l represents the photon flux measured by the 
ith detector. The y/s are assumed independent and that &,-,r, and {ay} are known 
non-negative constants. \i is also assumed to be independent of energy. 

To find a statistical estimate for the attenuation coefficient vector \x 
10 that is anatomically reasonable, a likelihood-based estimation approach is used. 
This is a natural choice since the likelihood is based on the statistical properties of 
the problem. The maximum likelihood (ML) approach also has good theoretical 
properties. It is asymptotically consistent, unbiased and efficient. The Poisson log 
likelihood for independent measurements is given by: 

15 L( M ) = t fclogCV- 1 ** + r,)- (V Iirt + 1})} (10) 

ignoring constant terms. 
Re gularization 

Without regularization, the Maximum Likelihood (ML) algorithm 
leads to noisy reconstruction. To reduce noise, it is possible to stop the algorithm 
20 when the reconstruction is visually appealing. Another approach is to pre-filter the 
data or post-filter the reconstruction. 

Regularization (penalized-likelihood) approaches for noise reduction 
have two important advantages. First, the penalty function improves the 
conditioning of the problem. Second, one can choose penalty functions with certain 
25 desirable properties such as edge preservation. 
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A general form for the regularizing penalty is the following: 

where i|/'s are potential functions acting on the soft constraints C\i«0 and K is the 
number of such constraints. Generally, the potential functions are symmetric, 
5 convex, non-negative and differentiable. Non-convex penalties can be useful for 
preserving edges, but are more difficult to analyze. One can think of the penalty 
as imposing a degree of smoothness or as a Bayesian prior. Both views are 
practically equivalent. 

Most commonly, penalty function penalize differences in the 
10 neighborhood of any particular pixel. They can be expressed as: 



15 



j= I keNj 



(12) 



where the weights w jk are 1 for orthogonal pixels and -j= for diagonal pixels and 
Nj is the pixel neighborhood. 

x 2 

A common penalty is the quadratic function, y (x) = — . This 
penalty is effective for noise reduction and is analytically simple, but can cause 
significant blurring, especially at edges. To preserve the edges, one can use a 
penalty that is less penalizing of large differences. One used herein is the Huber 
penalty : 



2 ' 



x<6 



S 2 



(13) 



20 A plot of both penalties is shown in Figure 4. 
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Combining the ML objective function with a penalty gives a 
penalized-likelihood (PL) objective function: 



<t><ji)=L(M)-fiR(ji). (14) 

The first term on the right in Equation (14) forces the estimator to 
5 match the measured data. The second term imposes some degree of smoothness 
leading to visually appealing images. The scalar parameter p controls the tradeoff 
between the two terms (or, alternatively, between resolution and noise). 

The goal of the reconstruction technique becomes to maximize (14) 
subject to certain object constraints such as non-negativity: 

. argmax 

io f ^0*™ (15) 



Ordered Subsets Transmission Algorithms 



The ideal statistical algorithm converges to the solution quickly (in 
a few iterations) and monotonically . It easily incorporates prior knowledge and 
constraints, accepts any type of system matrix and is parallelizable. No practical 
15 algorithm fits this description, and one settles for some compromise between the 
conflicting requirements. 

Two closely related statistical methods for X-ray CT, namely, the 
Ordered Subsets Transmission Reconstruction (OSTR) and the Ordered Subsets 
Penalized Weighted Least Squares (OS-PWLS) algorithms are described herein. 
20 OSTR models the sinogram data using Poisson statistics. The randomness in the 
measurement is a result of the photon emission and detection processes, electronic 
noise in the current-integrating detectors, as well as background radiation and 
scatter. Using the Poisson model is important in the low SNR case, since incorrect 
modeling can lead to reconstruction bias. At high SNR, the noise can be adequately 
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modeled as additive and Gaussian. The additive Gaussian model leads to a least- 
squares method. 

Both methods are based on the optimization transfer principle, which 
is discussed first. Then, the Poisson transmission model and associated algorithm 
5 is developed, including the idea or ordered subsets as a way to accelerate 
convergence. Although one could develop the least-squares approach from a 
Gaussian likelihood, it is derived as a quadratic approximation to the Poisson 
likelihood. 



Optimization Transfer Principle 

10 The optimization transfer principle is a very useful and intuitive 

principle that underlies many iterative methods, including the ones described herein. 
De Pierro introduced it in the context of inverse problems (emission tomography, 
to be specific). 

Often in iterative techniques, the goal is to maximize some objective 
15 function <S>(0) with respect to its argument 0. The objective G>(0) can be difficult to 
maximize. One resorts to replacing O with a surrogate function 4>(0;0 (n) ) that is 
easier to maximize at each iteration. Figure 5 illustrates the idea in one-dimension. 
The full utility of optimization transfer comes into play when the dimension of 0 is 
large, such as in tomography. 

20 The process is repeated iteratively, using a new surrogate function 

at each iteration. If the surrogate is chose appropriately, then the maximizer of 0(0) 
can be found. Sufficient conditions that ensure that the surrogate leads to a 
monotonic algorithm are known. 

Paraboloidal surrogates are used because they are analytically simple, 
25 and can be easily maximized. One can also take advantage of the convexity of these 
surrogates to parallelize the algorithm. 
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Separable Paraboloidal Surrogates 

Instead of maximizing the log likelihood, one can minimize the 
negative log likelihood, which is written as: 

(16) 

5 where h t (l) = bfi*+r r Yf0gdbje*+rd. Direct minimization of (16) leads to noisy 
estimates, so the likelihood is penalized by adding a roughness penalty. The 
problem becomes to find an estimate £ such that: 



argmin 



where 



10 0(/y) = - L(ji)+ fiR(fi). (18) 

The parameter p controls the tradeoff between the data-fit and penalty terms, and 
R(n) imposes a degree of smoothness on the solution. 

When r, * 0, the log likelihood is not concave and is difficult to 
maximize. A monotonic algorithm is possible if the optimization transfer principle 
15 is applied with paraboloidal surrogates. Paraboloidal surrogates are used such that 
the iterates ff monotonically decrease For that to occur, the surrogate <{> must 
satisfy the following "monotonicity" condition: 

<&(//)- d>(// w )< *K//";/0> V// > 0. (19) 

The following conditions are sufficient to ensure (19): 
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<D( A ") 



dfij 



, y = i,...,p 



(20) 



forfi 2> 0. 



Attention is restricted to different iable surrogates. Paraboloidal 
surrogates are convenient to use because they are easy to minimize. In (16), the 
5 likelihood function hfi) is replaced with the following quadratic surrogate: 



(21) 



where h x is the first derivative of h { . To ensure monotonicity, curvature q must be 
such that the surrogate satisfies the monotonicity conditions (19). Alternatively, to 
reduce computation, one can pre-compute the c/s prior to iterating. 

10 Replacing h x {[Aix]^ with q^fA/iJrJA^Ji) in (16) above will lead to a 

paraboloidal surrogate for the log likelihood. One can further take advantage of the 
nature of the surrogate to obtain a separable algorithm, i.e., an algorithm where all 
pixels are updated simultaneously. 

Rewrite the line integral: 



15 



(22) 



where 



£ a 0 = 1, V/,a tf > 0. 



7-1 



(23) 



Using the convexity of and writing /," = [Afi"\ , one gets: 
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3l 

\ a ij 



The overall separable surrogate for the log likelihood becomes: 



(24) 



(25) 



(26) 



where 



This formulation decouples the pixels. Each pixel effectively has its own cost 
function q y The q/s can be minimized for all pixels simultaneously, resulting in a 
parallelizable algorithm. 

10 A similar development can be pursued for the penalty term R(\l). 

Taking advantage of the convexity of the potential ftinction yields a separable 
penalty surrogate S(n;n n ). One now seeks to minimize the new separable global 
surrogate: 



(27) 



15 Since the surrogate is a separable paraboloid, it can be easily minimized by zeroing 
the first derivative. ( This leads to the following simultaneous update algorithm: 
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n 



j- l,-,p 



(28) 



where jl is some estimate of n, usually taken to be the current iterate \i n . The [•]+ 
operator enforces the non-negativity constant. From (21) and (25) the first and 
second derivatives of the surrogate are obtained: 



N 



dS 



1=1 



fi'M 



(29) 



dn) 



(30) 



10 



15 



q is the surrogate curvature and {a y } satisfies (23). To make the denominator in 
(28) small (and hence the step size large), one wants {a^ to be large. One also 
wants {a h ) to facilitate convergence, and to be independent of the current iterate so 
that it can be pre-computed. One convenient choice is: 



a„ = 



(31) 



It is possible that better choices exist. For the surrogate curvature 
c h one can use the optimal one derived in (24). The optimal curvature is iteration 
dependent. To save computations, one can use the following pre-computed 
approximation for the curvature: 



C: = h t log 

V 



(32) 



The pre-computed curvature may violate the conditions of 
monotonicity. It does, however, give an almost-monotonic algorithm, where the 
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surrogate becomes a quadratic approximation of the log likelihood. The pre- 
computed curvature seems to work well in practice, and the computational savings 
seem well worth the sacrifice. Major computational savings, however, come from 
the use of ordered subsets, discussed hereinbelow. 

5 Ordered Subsets 

Ordered subsets are usefiil when an algorithm involves a summation 
over sinogram indices (i.e. , a back projection). The basic idea is to break the set 
of sinogram angles into subsets, each of which subsamples the sinogram in the 
angular domain. The back projection process over the complete sinogram is 
10 replaced with successive back projections over the subsets of the sinogram. One 
iteration is complete after going through all of the subsets. 

Ordered subsets have been applied to emission tomography with a 
good degree of success. Improvements in convergence rate by a factor 
approximately equal to the number of subsets have been reported. Ordered subsets 

15 have also been used with transmission data for attenuation map estimation in 
SPECT. Ordered subsets where applied to the convex algorithm, and an increase 
in noise level with number of subsets have been reported. Ordered subsets have 
been used with the transmission EM algorithm and a cone-beam geometry. The 
OSTR algorithm was originally developed for attenuation correction in PET scans 

20 with considerable success. 

^ ' The choice of ordering the subsets is somewhat arbitrary, but it is 

preferably to order them "orthogonally." In such an arrangement, the projection 
corresponding to angles with maximum angular distance from previously chosen 
angles are used at each step. 

25 The cost of accelerating convergence with ordered subsets is loss of 

monotonicity. Hence, the term "convergence" is used loosely to mean that one gets 
visually acceptable reconstruction. Practically speaking, this loss of monotonicity 
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seems to make very little difference for the end result, especially if the algorithm 
is initialized with a reasonable starting image, such as an FBP reconstruction. 

Understanding the convergence properties of ordered subsets can 
provide insight into additional ways to accelerate convergence that may not be 
readily apparent. 

■ 

To summarize, the algorithm thus far flows as follows: 



10 



for each iteration n = 1 , . . . , niter 
- for each subset S - 1, .... M 

* compute h,([Ap]) 

* compute c,([Ap]) (or use pre-computed value) 
*L j =M£ ie! fl i jti l J=l,...,p 

* dj=M£ itS a l jy l c„ j=l,.:.,p 



15 



* = 



dj + fi 



d 2 S 



- end 
end 



20 



Scaling the denominator and numerator by the number of subsets 
ensures that the regularization parameter p remains independent of the number of 
subsets. This algorithm is known as the Ordered Subsets Transmission 
Reconstruction (OSTR) algorithm. 

OSTR combines the accuracy of statistical reconstruction with the 
accelerated rate of convergence that one gets from ordered subsets. The separability 
of the surrogates makes the algorithm easily parallelizable. The algorithm also 
naturally enforces the non-negativity constraint. The monotonicity property has 
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been satisfied, but that seems to hardly make a difference in practice if a reasonable 
starting image is used. 

Penalized Weighted Least Squares with Ordered Subsets 



5 high SNR scans, the Gaussian model is a reasonable approximation to the Poisson 
distribution. The Gaussian model leads to a simpler quadratic objection function 
and weighted-least-squares minimization. With high counts, PWLS leads to 
negligible bias and the simpler objective function reduces computation time. Figure 
6 illustrates how the quadratic approximation to the likelihood improves with count 
10 number. 



to the Poisson likelihood, which leads to a simpler objective function. The 
regularization term and the use of ordered subsets are retained. This variation of 
the method of the invention is called Ordered Subset Penalized Weight Least 
15 Squares (OS-PWLS). 



OSTR uses Poisson statistics to model the detection process. For 



The algorithm is reformulated by deriving a quadratic approximation 



For convenience, the negative log likelihood for transmission data is 



rewritten: 



N 




(33) 




(34) 



A 



20 



Taylor's expansion is applied to hflj) around some value / , and first 
and second order terms only are retained. 




2 



(35) 
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where h t and h t are the first and second derivatives of h ( . Assuming Y t > r„ one can 
estimate the line integral with: 



/, = log 



\Y-r. 



J 



(36) 



Substituting this estimate in (35) gives the following approximation for h t : 



(37) 



The first term in (37) is independent of /, and can be dropped. The 

(Y-r) 2 

weight is w, = . The new objective function is: 

Yf 



N 



(38) 



10 The subscript q indicates that this objective function is based on a quadratic 
approximation to the log likelihood. Subsequently, the subscript is dropped. The 
penalty term is also added. Minimizing this objective function over \i * 0 will lead 
to an estimator with negligible bias, since the number of counts is large. 



A separable surrogate for the objective function is almost 
15 immediately available. The terms inside the objective function summation are all 
convex. This convexity "trick" is explored yet one more time. Along the lines of 
(24) and (25), the surrogate for the PWLS objective function is: 



N p 

e,(^ n )=ZZ« 



' 2 



ft 
\ a ij 



(Mj-ftfl+Wl-l, 



(39) 



= Li j ■ 



(40) 
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The subscript q again emphasizes that this surrogate resulted from the quadratic 
approximation to the likelihood. It is dropped to simplify notation. A surrogate for 
the penalty function can be derived in a similar manner and added to (40). 

To use the iterative update (28) to minimize the surrogate, its first 
5 and second derivatives are computed, and the computational savings of the PWLS 
algorithm is demonstrated. 



= 2> # w f ([vw,) 



(41) 



dix) 

Unlike h. term in (29), the numerator (first derivative) in (41) 
10 involves no exponential terms and the denominator (second derivative) in (42) can 
be prercomputed and stored. The sum over sinogram indices can also be broken 
into sums over ordered subsets, further accelerating the algorithm. 

The iterative update is rewritten with the changes resulting from OS- 

PWLS: 



=z 



n u o vy i 



11= ll' 



f-1 



a 



v 



(42) 



dS 



15 



Ml ft 

a., 



d 2 S 



(43) 



Both OSTR and OS-PWLS have been described above for the 
monoenergetic case. The more realistic case of multi-energetic X-ray beam is 
described hereinbelow and a correction scheme is incorporated into OS-PWLS for 
the artifacts caused by the broad beam spectrum. 
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Statistical Reconstruction for P nlvenergetic X-ray CT 

One of the strengths of statistical methods is their applicability to 
different system and physical models. The CT transmission model is generalized 
hereinbelow to account for the broad spectrum of the beam. From the model, a cost 
5 function is derived and an iterative algorithm is developed for finding the minimizer 
of the cost function. 

Polvenergetic Statistical Model for CT 

A model for X-ray CT is described now that incorporates the energy 
dependence of the attenuation coefficient. A prior art algorithm could be applied 

10 to an image reconstructed with OS-PWLS. Instead, beam hardening correction is 
developed as an integral element of the statistical algorithm. The object is assumed 
to be comprised of K non-overlapping tissue types (this assumption may be 
generalized to allow for mixtures of tissues). For example, when K=2, one can use 
soft tissue and bone tissue classes, and when K=3, one can use soft tissue, bone, 

15 and a contrast agent, such as iodine. An iterative algorithm that generalizes OS- 
PWLS naturally emerges from the model. 

Assume that the if tissue classes are determined by pre-processing the 
data with soft tissue correction and then segmenting an initial reconstructed image. 
The object model of (2) is restricted to two spatial dimensions: 

20 M(x,y,E)=j^m k (E)p k (x >y y(x 9 y). (44) 



The system matrix is denoted with A = {a y } and the following definitions are made: 

R k = {set of pixels classified as tissue type k}> (45) 
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s}{p>[p\x,y)r k {x 9 y)dl, (47) 

v,{p)= sf ,...,*?). (48) 

The mass attenuation coefficients {m t (£)}£., of each of the K tissue types are 
assumed to be known. Discretization aside, from (1) the mean of the measured data 
5 along path L, is: 

= \l.(E)e- 1 " n * {E) ' :iP) dE 
J (49) 

= \l l {E)e' S! ' (E) -' ( - p) dE 

where = [m 1 (E),...,m k (E)] and the ' stands for vector transpose. The 
measurements are expressed as a function of the column vector quantity v, which has 
10 its elements the line integrals of the K different tissue densities. From knowledge 
of the X-ray spectrum, the values of ^(v,)and its gradient V Iftv,) are tabulated. 
In the discrete domain, 

*;(/>)=!>,/>,■ < 50 > 

The goal of the algorithm is to estimate the density coefficient vector 
15 />=/>„... pj. Rather than estimating K vector quantities of length />, each 
representing the density of one kind of tissue, the non-overlapping assumption of 
the tissue types enables one to keep the number of unknowns equal to as is the 
case in the monoenergetic model. This is possible when prior segmentation of the 
object is available. This can be obtained from a good FBP image, for example. 
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Polvenerpetic Model Cost Function 

m 

The Poisson log likelihood is set up in terms of the density p and the 
vector v,. To get a quadratic cost function, one follows a similar procedure to that 
described hereinabove, using the second-order Taylor's expansion. 

5 The function f ; represents the expected value of the measurement 

Y, at the ith detector. Using Y, in (10) gives the following negative log likelihood: 



-UP) = 2>,(v,(/0) ( 51 > 



Mv,(/>))= - ?logK(v,G>))+ r,J+ (I?(V/(P))+ '/) < 52 > 



The problem now is to find an estimate p such that: 

argmin 

10 /?= n (53) 



/>> 0 



where 



<b (j>) = - L(p) + fiRQ>) . (54) 

The regularization term can be treated exactly as before or it can be 
modified to avoid smoothing between different tissue types. For now, one focuses 
15 on the likelihood term and set the background term r, = 0. 

Suppose one can determine some initial estimate of v,fc), denoted 
v, = (3},...,s, K ). One expands h,(vj(p)) in a second-order Taylor series around v,: 



h,(Y)» MD+Vh^yv- ^(JJ- v,)'VV»,(v,)(v- v,). (55) 
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Taking the first and second derivatives of /i,.(v) = - Y, logi;(v) + Y,(v) , one gets 
the following for the first and second order terms of the Taylor expansion: 



V/i,(v) = 



1--= 



Y,(v). 



V^(v) 



(56) 



15 



V%(v) = 



1- * 



I Y,(v)) 



v 2 i;(v)+-^vi;(v)vi;(v). (57) 



5 The gradient VA, is a row vector and the Laplacian operator V z gives 

a Kx K matrix of partial derivatives. 

To simplify the algorithm and maintain the desirable property of 
separability, Y, is assumed to be close enough to l^(v ( .)for one to drop the first 
term on the right of (57). This also ensures that the resulting Hessian 
10 approximation is non-negative definite. 

In the Taylor expansion (55), the first term is constant and does not 
affect minimization, so it is dropped. The following is a quadratic approximation 
to the negative log likelihood: 



1 _ 2 (59) 

= .I 1 VA,(v / )(v/(p)-v / ) + — (V^g/XS/^)-*/)) • 

/= I 2.1; 



Substituting (50) into (59) and expanding the vector inner product yields: 
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121 



= 1 



£va(v,x*/o»)-^)+ 



21J 



(60) 



(61) 



where v k denotes the Jfcth element of the gradient vector. To simplify the above 
equation, the following definitions are made: 



10 



a: 



b.. A 
y = 



K 



Jfc=l 



A = {a,-,} is the geometrical system matrix. The matrix B = {b 0 } is 
a weighted system matrix, with the weights expressed as the non-zero elements of 
a diagonal matrix D(), to the left of A. The term Z t combines constants 
independent of /?. With the above definitions, expressing the line integrals explicitly 
in terms of the image pixels yields the following form of the cost function: 



N K 



t „ \ 



* „(p) ■ I S v Mi,) LVj(Pj - Pj) 

V 7=1 



2^ 



+ ftR(p). (62) 



This cost function is convex, so a separable surrogate and an iterative 
update are easily derived as described above. The results of the algorithm 
derivation are described below. 



15 



The separable paraboloidal surrogate for $ q (p) is given by: 
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p N K 



N K 



Q(p;p n ) = Z Z Z ? MKrfPj - Z Z v A(v,M/yV, 

j«l /=! k*\ 



+ Z Z ^7 a < 



(63) 



where 



Setting the point of linearization of the Taylor series at and 
5 evaluating the first and second derivative of Q at the same point gives: 



dp } 



N K 



= Z Zv/M(v,)+/? 



dS 



p-p 



(64) 



?Q(p;p") 



dp) 



A 1 ti, 



p=p 



/»i Yi <* t j fyj 



(65) 



10 



The overall algorithm is: 
initialize with p. 

for each iteration n = 1 , . . . , niter 
- for each subset 5 = 1,...,M 



compute 5* = ]£ a^pj for k= l,...,tf,v, = 



compute l^(v.), its gradient vector V Y,(v f ) and 
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compute byY v JXii) a u r j 
k~\ 

Zl 

K 

compute Nj = X £ v/ v **kC2 # ) 



compute 



Pj 



P=P 



d 2 S 



dp] 



p=p J 



- end 
end 



(66) 



10 



This algorithm globally converges to the minimizer of the cost 
function & q (p) when one subset is used, provided the penalty is chosen so that 4> q is 
strictly convex. When two or more subsets are used, it is expected to be monotone 
in the initial iterations. 



Referring now to Figure 10, there is illustrated in block diagram flow 
chart form, the method of the present invention. 

Initially, the spectrum of the incident X-ray beam is calibrated. 



15 



Then, the initial CT image is obtained. 



Segmentation of the CT image is then performed. 



Reprojection of the segmented image to calculate each component 
thickness is performed. 
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Then, measurement means and gradients are calculated as described 

above. 

Then, a cost function gradient is computed using back projection to 
yield the correction image as also described above. 

5 The correction image is then subtracted and non-negativity is 

enforced. 

The number of iterations is checked against a predetermined number 
or other criteria and if the iterative part of the method is complete, then the final 
corrected image is displayed. If not done, the iterative part of the method is re- 
10 entered at the reprojection step. At least some of the results obtained after 
subtraction may be used in the segmentation step as described herein as indicated 
by the dashed line from the "DONE" block. 

Realistically, many pixels are not comprised of one tissue only, but 
can be a mixture of several substances. This fact and the errors it causes in CT 

15 reconstruction, known as the partial volume effect, must be addressed for more 
accurate CT images. The method and apparatus described above are generalizable 
to the case of voxel mixtures by using fractional values in equation (46). It is 
possible to use histogram information to determine the value of the attenuation 
coefficient as a probabilistically weighted sum of several tissue types. Combining 

20 multi-energy measurements with tissue characteristics may also lead to more 
accurate "mixei models", where a pixel contains a tissue mixture. 

Al gorithm 

So far, beam hardening correction in the method of the invention 
depends on the availability of accurate classification of the different substances in 
25 the object. In simulated phantom studies, the bone/tissue distribution was known 
exactly. In a more realistic setting, such a classification would be available from 
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segmenting an initial image reconstructed with FBP. Using this segmentation map 
for all iterations may adversely affect the accuracy of the reconstruction. 

One promising alternative approach is to use joint likelihoods and 
penalties to estimate both pixel density values and tissue classes. In such an 

5 approach, the tissue classes are treated as random variables with a Markov random 
field model and are estimated jointly with the attenuation map. The joint likelihood 
will be a function of both the pixel density value and the pixel class. The joint 
penalties involve two parameters that balance the tradeoff between data fit and 
smoothness. In addition, joint penalties must account for the fact that pixels tend 

10 to have similar attenuation map values if the underlying tissue classes are the same, 
and vice versa. Such penalties would encourage smoothness in the same region but 
allow discontinuities between regions of different tissues. 

Com putation Time 

The long computation times of statistical methods hinder their use in 

* 

15 clinical X-ray CT applications. Significant acceleration by using ordered subsets 
have been demonstrated. 

From the algorithm design perspective, minor modifications may be 
made to accelerate convergence. Another possibility is to use a hybrid class of 
methods, which combines the faster early convergence rate of gradient methods with 
20 the faster ultimate linear convergence rate of steepest descent. 

> 

The most computationally expensive components of the iterative 
algorithm are back projection and forward projection, and there are algorithms that 
claim to perform these operations very quickly. It is possible that customized 
hardware may be used to perform the projections. Some recent work showed that 
25 readily available 2-D texture mapping hardware speeds up the Simultaneous 
Algebraic Reconstruction Technique (SART) to almost real-time realizations. 
SART involves forward and back projections much like OS-PWLS does. 
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Computation time may be reduced by using Fourier domain methods, 
where the 2-D Fourier transform of the image is assembled from its projections. 
The practicality of this approach depends on the availability of fast gridding 
methods in the Fourier domain. 

5 Scatter 

Scatter is a major problem in cone-beam CT, where it can range from 
40% up to 200% of the unscattered data. Collimation reduces scatter, but 
collimating flat-panel detectors is challenging. Among several factors that affect the 
performance of a cone-beam, flat-panel detector computed tomography system, 
10 scatter was shown to degrade the detector quantum efficiency (DQE) and to 
influence the optimal magnification factor. Larger air gaps were needed to cope 
with high scatter, especially if imaging a large FOV. 

There are generally two approaches to dealing with the problem of 
scatter. Scatter can either be physically removed (or reduced) before detection or 
15 can be numerically estimated and its effect compensated for. Ways to physically 
remove scatter include air gaps and grids, but are not very practical once flat-panel 
detectors are used, due to the small size of detector pixels. 

There are several approaches to numerically estimate scatter. For 
example, idle detectors (those outside FOV) can be used to provide a scatter 
20 estimate. Such measurements can be combined with analytic models for scatter that 
depend, among other factors, on the energy of the radiation used, the volume of 
scattering material and system geometry. 

There have also been encouraging attempts at incorporating scatter 
estimation/correction into statistical reconstruction. For instance, ordered subsets 
25 EM methods with scatter correction gave superior results to FBP with scatter 
correction in emission CT. In another work, the maximum likelihood EM 
algorithm successfully improved contrast and SNR of digital radiology images by 
incorporating a convolution-based scatter estimate. This is particularly relevant 



-36- 



WO 02/067201 



PCT/US01/04894 



since cone-beam transmission CT with flat-panel detectors is in many ways similar 
to digital radiography. It is, in essence, a rotating digital radiography system. 

A scatter estimate may be incorporated in the model of the CT 
problem as well as in the reconstruction algorithm using the ri terms above. One 
5 approach to scatter estimation and correction is a numerical one, i.e., ways to 
physically eliminate scatter will not be considered. This makes the approach 
portable to different systems, and less costly. 

System Desien 

The statistical measurement model has the potential to be extended 
10 to take into account the time dimension when imaging the heart, and thus free the 
designer from synchronization constraints. Moreover, statistical reconstruction 
eliminates the need for rebinning and interpolation. This may lead to higher helical 
scanning pitch and cardiac imaging with good temporal and axial resolutions. 

Conclusion 

15 In conclusion, a framework has been described above for using 

statistical methods to reconstruct X-ray CT images from poly energetic sources. 

While the best mode for carrying out the invention has been 
described in detail, those familiar with the art to which this invention relates will 
recognize various alternative designs and embodiments for practicing the invention 
20 as defined by the following claims. 
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WHAT IS CLAIMED IS: 



1 . A method for statistically reconstructing a polyenergetic X-ray 
computed tomography image to obtain a corrected image, the method comprising: 

providing a computed tomography initial image produced by a single 
5 X-ray CT scan having a polyenergetic source spectrum wherein the initial image has 
components of materials which cause beam hardening artifacts; 

separating the initial image into different sections to obtain a 
segmented image; and 

calculating a series of intermediate corrected images based on the 
10 segmented image utilizing a statistical algorithm which accounts for the 
polyenergetic source spectrum and which converges to obtain a final corrected 
image which has significantly reduced beam hardening artifacts. 

2. The method as claimed in claim 1 wherein the step of 
calculating includes the steps of calculating a gradient of a cost function having an 

15 argument and utilizing the gradient to minimize the cost function with respect to its 
argument. 



3. The method as claimed in claim 2 wherein the step of 
calculating the gradient includes the step of back projecting. 

4. The method as claimed in claim 2 wherein the step of 
20 calculating the gradient includes the step of calculating thicknesses of the 

components . 

5. The method as claimed in claim 4 wherein the step of 
calculating thicknesses includes the step of reprojecting the segmented image. 

6. The method as claimed in claim 4 wherein the step of 
25 calculating the gradient includes the step of calculating means of data along paths 

and gradients based on the thicknesses of the components. 



-38- 



WO 02/067201 PCT/US01/04894 



7. The method as claimed in claim 2 wherein the argument is 
density of the materials at each image voxel. 

8. The method as claimed in claim 1 further comprising 
calibrating the spectrum of the X-ray CT scan. 

5 9. The method as claimed in claim 1 further comprising 

displaying the final corrected image. 

10. The method as claimed in claim 2 wherein the step of 
calculating the gradient includes the step of utilizing ordered subsets to accelerate 
convergence of the algorithm. 

10 1 1 . The method as claimed in claim 2 wherein the cost function 

has a regularizing penalty term. 

12. An image reconstructor apparatus for statistically 
reconstructing a polyenergetic X-ray computed tomography image to obtain a 
corrected image, the apparatus comprising: 

15 means for providing a computed tomography initial image produced 

by a single X-ray CT scan having a polyenergetic source spectrum wherein the 
initial image has components of materials which cause beam hardening artifacts; 

means for separating the initial image into different sections to obtain 

a segmented image; and 
20 means for calculating a series of intermediate corrected images based 

on the segmented image utilizing a statistical algorithm which accounts for the 
polyenergetic source spectrum and which converges to obtain a final corrected 
image which has significantly reduced beam hardening artifacts. 

13. The apparatus as claimed in claim 12 wherein the means for 
25 calculating includes means for calculating a gradient of a cost function having an 

argument and means for utilizing the gradient to minimize the cost function with 
respect to its argument. 
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14. The apparatus as claimed in claim 13 wherein the means for 
calculating the gradient includes means for back projecting. 

15. The apparatus as claimed in claim 13 wherein the means for 
calculating the gradient includes means for calculating thicknesses of the 

5 components. 

16. The apparatus as claimed in claim 15 wherein the means for 
calculating thicknesses includes means for reprojecting the segmented image. 

17. The apparatus as claimed in claim 15 wherein the means for 
calculating the gradient includes means for calculating means of data along paths 

10 and gradients based on the thicknesses of the components. 

18. The apparatus as claimed in claim 13 wherein the argument 
is density of the materials at each image voxel. 

19. The apparatus as claimed in claim 12 wherein the means for 
calculating the gradient includes means for utilizing ordered subsets to accelerate 

15 convergence of the algorithm. 

20. The apparatus as claimed in claim 13 wherein the cost 
function has a regularizing penalty term. 
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